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Abstract. Fluctuations in the abundance of molecules in the living cell may 
affect its growth and well being. For regulatory molecules (e.g., signaling proteins 
or transcription factors), fluctuations in their expression can affect the levels of 
downstream targets in a network. Here, we develop an analytic framework to 
investigate the phenomenon of noise correlation in molecular networks. Specifically, 
we focus on the metabolic network, which is highly inter-linked, and noise properties 
may constrain its structure and function. Motivated by the analogy between the 
dynamics of a linear metabolic pathway and that of the exactly soluable linear 
queueing network or, alternatively, a mass transfer system, we derive a plethora 
of results concerning fluctuations in the abundance of intermediate metabolites in 
various common motifs of the metabolic network. For all but one case examined, we 
find the steady-state fluctuation in different nodes of the pathways to be effectively 
uncorrelated. Consequently, fluctuations in enzyme levels only affect local properties 
and do not propagate elsewhere into metabolic networks, and intermediate metabolites 
can be freely shared by different reactions. Our approach may be applicable to study 
metabolic networks with more complex topologies, or protein signaling networks which 
are governed by similar biochemical reactions. Possible implications for bioinformatic 
analysis of metabolimic data are discussed. 
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Due to the limited number of molecules for typical molecular species in microbial 
cells, random fluctuations in molecular networks are common place and may play 
important roles in vital cellular processes. For example, noise in sensory signals 
can result in pattern formation and collective dynamics [I], and noise in signaling 
pathways can lead to cell-to-cell variability [2J. Also, stochasticity in gene expression 
has implications on cellular regulation OH] and may lead to phenotypic diversity [3 [6], 
while fluctuations in the levels of (toxic) metabolic intermediates may reduce metabolic 
efficiency [7j and impede cell growth. 

In the past several years, a great deal of experimental and theoretical efforts have 
focused on the stochastic expression of individual genes, at both the translational and 
transcriptional levels [8], [9J [10] . The effect of stochasticity on networks has been studied 
in the context of small, ultra-sensitivie genetic circuits, where noise at a circuit node (i.e., 
a gene) was shown to either attenuate or amplify output noise in the steady state [TT|[T2"]. 
This phenomenon — termed 'noise propagation' — make the steady-state fluctuations 
at one node of a gene network dependent in a complex manner on fluctuations at other 
nodes, making it difficult for the cell to control the noisiness of individual genes of 
interest p3]. Several key questions which arise from these studies of genetic noise 
include (i) whether stochastic gene expression could further propagate into signaling 
and metabolic networks through fluctuations in the levels of key proteins controlling 
those circuits, and (ii) whether noise propagation occurs also in those circuits. 

Recently, a number of approximate analytical methods have been applied to 
analyze small genetic and signaling circuits; these include the independent noise 
approximation [Til [T5| [16] . the linear noise approximation [Tll[TT] . and the self-consistent 
field approximation [18]. Due perhaps to the different approximation schemes used, 
conflicting conclusions have been obtained regarding the extent of noise propagation 
in various networks (see, e.g., [1TJ . ) Moreover, it is difficult to extend these studies to 
investigate the dependences of noise correlations on network properties, e.g., circuit 
topology, nature of feedback, catalytic properties of the nodes, and the parameter 
dependences (e.g., the phase diagram). It is of course also difficult to elucidate these 
dependences using numerical simulations alone, due to the very large degrees of freedoms 
involved for a network with even a modest number of nodes and links. 

In this study, we describe an analytic approach to characterize the probability 
distribution for all nodes of a class of molecular networks in the steady state. Specifically, 
we apply the method to analyze fluctuations and their correlations in metabolite 
concentrations for various core motifs of the metabolic network. The metabolic network 
consists of nodes which are the metabolites, linked to each other by enzymatic reactions 
that convert one metabolite to another. The predominant motif in the metabolic 
network is a linear array of nodes linked in a given direction (the directed pathway), 
which are connected to each other via converging pathways and diverging branch points 
|19j . The activities of the key enzymes are regulated allosterically by metabolites 
from other parts of the network, while the levels of many enzymes are controlled 
transcriptionally and are hence subject to deterministic as well as stochastic variations 
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in their expressions [20J . To understand the control of metabolic network, it is important 
to know how changes in one node of the network affect properties elsewhere. 

Applying our analysis to directed linear metabolic pathways, we predict that 
the distribution of molecule number of the metabolites at intermediate nodes to 
be statistically independent in the steady state, i.e., the noise does not propagate. 
Moreover, given the properties of the enzymes in the pathway and the input flux, we 
provide a recipe which specifies the exact metabolite distribution function at each node. 
We then show that the method can be extended to linear pathways with reversible 
links, with feedback control, to cyclic and certain converging pathways, and even to 
pathways in which flux conservation is violated (e.g., when metabolites leak out of 
the cell). We find that in these cases correlations between nodes are negligable or 
vanish completely, although nontrivial fluctuation and correlation do dominate for a 
special type of converging pathways. Our results suggest that for vast parts of the 
metabolic network, different pathways can be coupled to each other without generating 
complex correlations, so that properties of one node (e.g., enzyme level) can be changed 
over a broad range without affecting behaviors at other nodes. We expect that the 
realization of this remarkable property will shape our understanding of the operation 
of the metabolic network, its control, as well as its evolution. For example, our results 
suggest that correlations between steady-state fluctuations in different metabolites bare 
no information on the network structure. In contrast, temporal propagation of the 
response to an external perturbation should capture - at least locally - the morphology 
of the network. Thus, the topology of the metabolic network should be studied during 
transient periods of relaxation towards a steady-state, and not at steady-state. 

Our method is motivated by the analogy between the dynamics of biochemical 
reactions in metabolic pathways and that of the exactly solvable queueing systems [16] 
or, alternatively, as mass transfer systems [221 IS]- Our approach may be applicable 
also to analyzing fluctuations in signaling networks, due to the close analogy between 
the molecular processes underlying the metabolic and signaling networks. To make our 
approach accessible to a broad class of circuit modelers and bioengineers who may not 
be familiar with nonequilibrium statistical mechanics, we will present in the main text 
only the mathematical results supported by stochastic simulations, and defer derivations 
and illustrative calculations to the Supporting Materials. While our analysis is general, 
all examples are taken from amino-acid biosynthesis pathways in E. Coli |24j . 

1. Individual Nodes 

1.1. A molecular Michaelis-Menton model 

In order to set up the grounds for analyzing a reaction pathway and to introduce our 
notation, we start by analyzing fluctuations in a single metabolic reaction catalyzed by 
an enzyme. 

Recent advances in experimental techniques have made it possible to track the 
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enzymatic turnover of a substrate to product at the single- molecule level [2SIEZ], and 
to study instantaneous metabolite concentration in the living cell [28]. To describe 
this fluctuation mathematically, we model the cell as a reaction vessel of volume V, 
containing m substrate molecules (S) and N E enzymes (E) . A single molecule of 5" can 
bind to a single enzyme E with rate k + per volume, and form a complex, SE. This 
complex, in turn, can unbind (at rate kS) or convert 5* into a product form, P, at rate 
k 2 . This set of reactions is summarized by 

S + E^SE^P + E . (1) 

fc_ 

Analyzing these reactions within a mass-action framework — keeping the substrate 
concentration fixed, and assuming fast equilibration between the substrate and the 
enzymes (k± ^> k 2 ) - - leads to the Michaelis-Menten (MM) relation between the 
macroscopic flux c and the substrate concentration [S] = m/V : 

c = v mBX [S\/([S\+K M ) , (2) 

where Km = k_/k + is the dissociation constant of the substrate and the enzyme, 
and f max = k 2 [E] is the maximal flux, with [E] = N E /V being the total enzyme 
concentration. 

Our main interest is in noise properties, resulting from the discreteness of molecules. 
We therefore need to track individual turnover events. These are described by the 
turnover rate w m , defined as the inverse of the mean waiting time per volume between 
the (uncorrelatecyj) synthesis of one product molecule to the next. Assuming again 
fast equilibration between the substrate and the enzymes, the probability of having 
N SE complexes given m substrate molecules and N E enzymes is simply given by the 
Boltzmann distribution, 

. K- N *>» m\N E \ 

p{N SE \m, N E ) = — —— — — r 3 

Z m ,N E N SE \(m - N SE )\{N E - N S E)i 

for N$ E < N E and m. Here K~ l = Vk + /k_ is the Boltzmann factor associated 
with the formation of an SE complex, and the Z m ,N E takes care of normalization (i.e., 
chosen such that ^2n se p(Ns E \m, N e ) = 1.) Under this condition, the turnover rate 
w m = y ^2 Ns E ■ p(Ns E \m, N E ) is given approximately by 

Tfl 

W " = V - rn + (K + N E -l) + ° (A " 3) < (4) 
with v max = k 2 N E /V; see Supp. Mat. We note that for a single enzyme (N E = 1), one 
has w m = v max m/(m + K), which was derived and verified experimentally [271 ESJ. 

1.2. Probability distribution of a single node 

In a metabolic pathway, the number of substrate molecules is not kept fixed; rather, 
these molecules are synthesized or imported from the environment, and at the same time 

% We note in passing that some correlations do exist - but not dominate - in the presence of "dynamical 
disorder" [27] , or if turnover is a multi-step process [29l [30] . 
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turned over into products. We consider the influx of substrate molecules to be a Poisson 
process with rate c. These molecules are turned into product molecules with rate w m 
given by Eq. (H]). The number of substrate molecules is now fluctuating, and one can ask 
what is the probability 7r(m) of finding m substrate molecules at the steady-state. This 
probability can be found by solving the steady-state Master equation for this process 
(see Supp. Mat.), yielding 



7r(m) 



where z = c/t> max [31]. The form of this distribution is plotted in supporting figure 1 
(solid black line). As expected, a steady state exists only when c < v max . Denoting the 
steady-state average by angular brackets, i.e., (x m ) = J2m x m 1T ( rn )^ the condition that 
the incoming flux equals the outgoing flux is written as 

C={Wm) = V - S + (K + N E ) > (6) 
where s = (m). 

Comparing this microscopically-derived flux-density relation with the MM relation 
([2]) using the obvious correspondence [S] = s/V, we see that the two are equivalent 
with Km = {K + Ne) /V. Note that this microscopically-derived form of MM constant 
is different by the amount [E] from the commonly used (but approximate form) 
Km = K/V, derived from mass-action. However, for typical metabolic reactions, 
Km ~ 10 — 1000 fiM [21] while [E] is not more than 1000 molecules in a bacterium 
cell (~ lfiM); so the numerical values of the two expressions may not be very different. 

We will characterize the variation of substrate concentration in the steady-state by 
the noise index 

n 2 = = t ' max (7) 
ls ~s 2 c-{K + N E )' K) 

where a 2 is the variance of the distribution n{m). Since c < t> max and increases with 
s towards 1 (see Eq. [6]), r] s decreases with the average occupancy s as expected. It 
is bound from below by 1/y/K + Ne, which can easily be several percent. Generally, 
large noise is obtained when the reaction is catalyzed by a samll number of high-affinity 
enzymes (i.e., for low K and Ne). 



2. Linear pathways 

2.1. Directed pathways 

We now turn to a directed metabolic pathway, where an incoming flux of substrate 
molecules is converted, through a series of enzymatic reactions, into a product flux [19] . 
Typically, such a pathway involves the order of 10 reactions, each takes as precursor the 
product of the preceding reaction, and frequently involves an additional side-reactant 
(such as a water molecule or ATP) that is abundant in the cell (and whose fluctuations 
can be neglected). As a concrete example, we show in figure [TJa) the tryptophan 
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Figure 1. Linear biosynthesis pathway, (a) Tryptophan biosynthesis pathway in E. 
Colt, (b) Model for a directed pathway. Dashed lines depict end-product inhibition. 
Abbreviations: CPAD5P, 1-O-Carboxyphenylamino l-deoxyribulose-5-phosphate; 
NPRAN, N-5-phosphoribosyl-anthranilate; IGP, Indole glycerol phosphate; PPL 
Pyrophosphate; PRPP, 5-Phosphoribosyl-l-pyrophosphate; T3P1, Glyceraldehyde 3- 
phosphate. 



biosynthesis pathway of E. Coli [21], where an incoming flux of chorismate is converted 
through 6 directed reactions into an outgoing flux of tryptophan, making use of several 
side-reactants. Our description of a linear pathway includes an incoming flux c of 
substrates of type Si along with a set of reactions that convert substrate type Si to 
S i+ i by enzyme E( (see figure [1(b)) with rate = v i m i /(m i + Ki — 1) according to 
Eq. (HI). We denote the number of molecules of intermediate Si by mj, with mi for the 
substrate and for the end-product. The superscript (i) indicates explicitly that the 
parameters Vi = /V and Ki = (K^ + Ng ) describing the enzymatic reaction 

Si — > Si + \ are expected to be different for different reactions. 

The steady-state of the pathway is fully described by the joint probability 
distribution 7r(mi,m2, . . . ,itll) of having rrii molecules of intermediate substrate type 
Si. Surprisingly, this steady-state distribution is given exactly by a product measure, 



where 7Tj(m) is as given in Eq. ([5]) (with K + Ne replaced by Ki and z by Zi = c/vi), 
as we show in Supp. Mat. This result indicates that in the steady state, the number of 
molecules of one intermediate is statistically independent of the number of molecules of 
any other substrat j§|. The result has been derived previously in the context of queueing 
networks |46] . and of mass-transport systems [17]. Either may serve as a useful analogy 
for a metabolic pathway. 

Since the different metabolites in a pathway are statistically decoupled in the steady 
state, the mean Sj = (mi) and the noise index rj^. = c~ x VijKi can be determined by 
Eq. (GO) individually for each node of the pathway. It is an interesting consequence 
of the decoupling property of this model that both the mean concentration of each 
substrate and the fluctuations depend only on the properties of the enzyme immediately 
downstream. While the steady-state flux c is a constant throughout the pathway, the 



L 




(8) 



§ We note, however, that short-time correlations between metabolites can still exist, and may be probed 
for example by measuring two-time cross-correlations; see discussion at the end of the text. 
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Figure 2. Noise in metabolite molecular number (r] a — cr s /s) for different pathways. 
Monte-Carlo simulations (bars) are compared with the analytic prediction (symbols) 
obtained by assuming decorrelation for different nodes of the pathways. The structure 
of each pathway is depicted under each panel. Parameter values were chosen randomly 
such that 10 3 < Ki < 10 4 and c < uj < 10c. Similar decorrelation was obtained for 
100 different random choices of parameters, and for 100 different sets with Ki 10-fold 
smaller (data not shown). The effect on the different metabolites of a change in the 
velocity of the first reaction, v\ — 1.1c (dark gray)-^ 5c (light gray), is demonstrated. 
Similar results are obtained for changes in K\ (data not shown.) (a) Directed pathway. 
Here the decorrelation property is exact, (b) Directed pathway with two reversible 



reactions. 

and -K^~4 



For these reactions, 



8.4, 6.9c; v^ A = 1.6,3.7,c;^ 4 
7700,3700. (c) Linear dilution of metabolites. Here /3/c ■ 



U 3A 



■■ 2500,8000 
1/100. (d) 

End-product inhibition,where the influx rate is given by a = cq [1 + (mjj/Ki)] 1 with 
Kj = 1000. (e) Diverging pathways. Here metabolite 4 is being processed by two 
enzymes (with different affinities, K 1 = 810, K 11 = 370) into metabolites 5 and 7, 
resp. (f) Converging pathways. Here two independent 3-reaction pathways , with 
fluxes c and c' = c/2, produce the same product, Si. 



parameters v j and Ki can be set separately for each reaction by the copy-number and 
kinetic properties of the enzymes (provided that c < Vi). Hence, for example, in a case 
where a specific intermediate may be toxic, tuning the enzyme properties may serve to 
decrease fluctuations in its concentration, at the price of a larger mean. To illustrate 
the decorrelation between different metabolites, we examine the response of steady-state 
fluctuations to a 5-fold increase in the enzyme level [Ex]. Typical time scale for changes 
in enzyme level much exceeds those of the enzymatic reactions. Hence, the enzyme level 
changes may be considered as quasi-steady state. In figure [2]^a) we plot the noise indices 
of the different metabolites. While noise in the first node is significantly reduced upon 
a 5-fold increase in [E\], fluctuations at the other nodes are not affected at all. 

2.2. Reversible reactions 

The simple form of the steady-state distribution (El) for the directed pathways may 
serve as a starting point to obtain additional results for metabolic networks with more 
elaborate features. We demonstrate such applications of the method by some examples 
below. In many pathways, some of the reactions are in fact reversible. Thus, a 
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metabolite Si may be converted to metabolite S i+ i with rate ^i ax ^i/(^i + Kf) or 
to substrate S^i with rate v^^rrii/ '(mi + K~). One can show — in a way similar to 
Ref. [32] — that the decoupling property ([5]) holds exactly only if the ratio of the two 
rates is a constant independent of m^, i.e. when Kf = K~ . In this case the steady state 
probability is still given by ([5]), with the local currents obeying 



This is nothing but the simple fact that the overall flux is the difference between the 
local current in the direction of the pathway and that in the opposite direction. 

In general, of course, Kf ^ K~ . However, we expect the distribution to be given 
approximately by the product measure in the following situations: (a) Kf ~ K~\ (b) 
the two reactions are in the zeroth-order regime, s ^> K i ; (c) the two reactions are in 

the linear regime, s -C Kf. In the latter case Eq. ([9]) is replaced by — -jr^-Zk+i = c . 

K i K i+i 

Taken together, it is only for a narrow region (i.e., Si ~ Ki) where the product measure 
may not be applicable. This prediction is tested numerically, again by comparing two 
pathways (now containing reversible reactions) with 5-fold difference in the level of 
the first enzyme. From figur42](b), we see again that the difference in noise indices 
exist only in the first node, and the computed value of the noise index at each node 
is in excellent agreement with predictions based on the product measure (symbols). 
Similar decorrelation was obtained for 100 different random choices of parameters, and 
for 100 different sets with Ki 10-fold smaller (data not shown). 

2.3. Dilution of intermediates 

In the description so far, we have ignored possible catabolism of intermediates or dilution 
due to growth. This makes the flux a conserved quantity throughout the pathway, and 
is the basis of the flux-balance analysis [32]. One can generalize our framework for the 
case where flux is not conserved, by allowing particles to be degraded with rate u m . 
Suppose, for example, that on top of the enzymatic reaction a substrate is subjected 
to an effective linear degradation, u m = (3m. This includes the effect of dilution due 
to growth, in which case (3 = ln(2)/(mean cell division time), and the effect of leakage 
out of the cell. As before, we first consider the dynamics at a single node, where the 
metabolite is randomly produced (or transported) at a rate cq. It is straightforward to 
generalize the Master equation for the microscopic process to include u m , and solve it 
in the same way. With w m as before, the steady state distribution of the substrate pool 
size is then found to be 



where (a) m = a(a + 1) ■ ■ ■ (a + m — 1). This form of 7r(m) allows one to easily calculate 
moments of the molecule number from the partition function Z as in equilibrium 
statistical mechanics, e.g. s = (m) = c dZ/dc , and thence the outgoing flux, c = c — /3s. 
Using the fact that Z can be written explicitly in terms of hypergeometric functions, 



V- Zi - v i+1 z i+1 = c . 



(9) 




(10) 
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we find that the noise index grows with (3 as rjl ~ v/(Kc ) + (3/c$. The distribution 
function is given in supporting figure 1 for several values of j3. 

Generalizing the above to a directed pathway, we allow for (3, as well as for t> max 
and K, to be 2-dependent. The decoupling property ([8]) does not generally hold in the 
non- conserving case [33]. However, in this case the stationary distribution still seems 
to be well approximated by a product of the single-metabolite functions 71* (m) of the 
form (1101). with c /j3 — > This is supported again by the excellent agreement 

between noise indices obtained by numerical simulations and analytic calculations using 
the product measure Ansatz, for linear pathways with dilution of intermediates; see 
figure Etc). In this case, change in the level of the first enzyme does "propagate" to 
the downstream nodes. But this is not a "noise propagation" effect, as the mean fluxes 
(cj) at the different nodes are already affected. (To illustrate the effect of leakage, the 
simulation used parameters that corresponded to a huge leakage current which is 20% 
of the flux. This is substantially larger than typical leakage encountered, say due to 
growth-mediated dilution, and we do not expect propagation effects due to leakage to 
be significant in practice.) 

3. Interacting pathways 

The metabolic network in a cell is composed of pathways of different topologies. While 
linear pathways are abundant, one can also find circular pathways (such as the TCA 
cycle), converging pathways and diverging ones. Many of these can be thought of as a 
composition of interacting linear pathways. Another layer of interaction is imposed 
on the system due to the allosteric regulation of enzyme activity by intermediate 
metabolites or end products. To what extent can our results for a linear pathway 
be applied to these more complex networks? Below we address this question for a few of 
the frequently encountered cases. To simplify the analysis, we will consider only directed 
pathways and suppress the dilution/leakage effect. 

3.1. Cyclic pathways 

We first address the cyclic pathway, in which the metabolite Sl is converted into Si 
by the enzyme El- Borrowing a celebrated result for queueing networks [31] and mass 
transfer models [35], we note that the decoupling property (EJ) described above for the 
linear directed pathway also holds exactly even for the cyclic pathways^. This result is 
surprising mainly because the Poissonian nature of the "incoming" flux assumed in the 
analysis so far is lost, replaced in this case by a complex expression, e.g., uiml • jtlW. 
In an isolated cycle the total concentration of the metabolites, s tot - and not the 

|| In fact, the decoupling property holds for a general network of directed single-substrate reactions, 
even if the network contains cycles. 
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flux - is predetermined. In this case, the flux c is give by the solution to the equation 



Note that this equation can always be satisfied by some positive c that is smaller than 
all Vi's. In a cycle that is coupled to other branches of the network, flux may be governed 
by metabolites going into the cycle or taken from it. In this case, flux balance analysis 
will enable determination of the variables which specify the probability distribution 



3.2. End-product inhibition 

Many biosynthesis pathways couple between supply and demand by a negative feedback 
[2"H [T5] , where the end-product inhibits the first reaction in the pathway or the transport 
of its precursor; see, e.g., the dashed lines in figure [H In this way, flux is reduced when 
the end-product builds up. In branched pathways this may be done by regulating 
an enzyme immediately downstream from the branch-point, directing some of the flux 
towards another pathway. 

To study the effect of end-product inhibition, we consider inhibition of the inflow 
into the pathway. Specifically, we model the probability at which substrate molecules 
arrive at the pathway by a stochastic process with exponentially-distributed waiting 
time, characterized by the rate a{mi) = Cq [l + (m^/i^/) 71 ] , where Cq is the maximal 
influx (determined by availability of the substrate either in the medium or in the 
cytoplasm), trl is the number of molecules of the end-product (Sl), Kj is the 
dissociation constant of the interaction between the first enzyme E and Sl, and h is a 
Hill coefficient describing the cooperativity of interaction between E and Sl- Because 
mi is a stochastic variable itself, the incoming flux is described by a nontrivial stochastic 
process which is manifestly non-Poissonian. 

The steady-state flux is now 



This is an implicit equation for the flux c, which also appears in the right-hand side of 
the equation through the distribution 7r(mi, ...,mi). 

By drawing an analogy between feedback-regulated pathway and a cyclic pathway, 
we conjecture that metabolites in the former should be effectively uncorrelated. The 
quality of this approximation is expected to become better in cases where the ration 
between the influx rate a{m L ) and the outflux rate w mL is typically idependent. 
Under this assumption, we approximate the distribution function by the product 
measure ([8]), with the form of the single node distributions given by ([5]). Note that the 
conserved flux then depends on the properties of the enzyme processing the last reaction, 
and in general should be influenced by the fluctuations in the controlling metabolite. In 
this sense, these fluctuations propagate throughout the pathway at the level of the mean 




(11) 



(El). 




(12) 
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flux. This should be expected from any node characterized by a high control coefficient 

Using this approximate form, Eq. ( 1121 ) can be solved self-consistently to yield c(cq), 
as is shown explicitly in Supp. Mat. for h = 1. The solution obtained is found to be 
in excellent agreement with numerical simulation (Supporting figure 2a). The quality 
of the product measure approximation is further scrutinized by comparing the noise 
index of each node upon increasing the enzyme level of the first node 5-fold. Figure [2]^c) 
shows clearly that the effect of changing enzyme level does not propagate to other nodes. 
While being able to accurately predict the flux and mean metabolite level at each node, 
the predictions based on the product measure are found to be under-estimating the 
noise index by up to 10% (compare bars and symbols). We conclude that in this case 
correlations between metabolites do exist, but not dominate. Thus analytic expressions 
dervied from the decorrelation assumption can be useful even in this case (see supporting 
figure 2b). 



3. 3. Diverging pathways 

Many metabolites serve as substrates for several different pathways. In such cases, 
different enzymes can bind to the substrate, each catabolizes a first raction in a different 
pathway. Within our scheme, this can be modeled by allowing for a metabolite Si to be 
converted to metabolite S\ with rate w^. = v I m i /(m i + K 1 — 1) or to metabolite Sf 1 
with rate = v IT mi/ '(mi + K 11 — 1). The paramters f 1 ' 11 and K I,JI characterize the 
two different enzymes. 

Similar to the case of reversible reactions, the steady-state distribution is given 
exactly by a product measure only if w^./w^. is a constant, independent of m; (namely 
when K 1 = K 11 ). Otherwise, we expect it to hold in a range of alternative scenarios, 
as described for reversible pathways. 

Considering a directed pathway with a single branch point, the distribution ([5]) 
describes exactly all nodes upstream of that point. At the branchpoint, one replaces 
w m by w m = + w mi to obtain the distribution function 

( ) ^ (K I )rn(K II )m , . 

{ ' Z ml^iW 1 + K 1I v 1I )/(v 1 +v 1I )) m ' 1 ' 

From this distribution one can obtain the fluxes going down each one of the two 
branching pathway, c 1,11 = Y2 w m Il7r ( m )- Both fluxes depend on the properties of 
both enzymes, as can be seen from (fl3l). and thus at the branch-point the two pathways 
influence each other [36]. Moreover, fluctuations at the branch point to propagate into 
the branching pathways already at the level of the mean flux. This is consistent with the 
fact that the branch node is expected to be characterized by a high control coefficient 

While different metabolite upstream and including the branch point are 
uncorrelated, this is not exactly true for metabolites of the two branches. Nevertheless, 
since these pathways are still directed, we further conjecture that metabolites in the two 
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Figure 3. Converging pathways, (a) Glycine is synthesized in two independent 
pathways, (b) Citrulline is synthesized from products of two pathways. Abbreviations: 
2A30, 2-Amino-3-oxobutanoate. 



branching pathways can still be described, independently, by the probability distribution 
([5]), with c given by the flux in the relevant branch, as calculated from (fT3l). Indeed, the 
numerical results of figure ^e) strongly support this conjecture. We find that changing 
the noise properties of a metabolite in the upstream pathway do not propagte to those 
of the branching pathways. 

3.4- Converging pathways - combined fluxes 

We next examine the case where two independent pathways result in synthesis of the 
same product, P. For example, the amino acid glycine is the product of two (very short) 
pathways, one using threonine and the other serine as precursors (figure [3]^a)) |24j . With 
only directed reactions, the different metabolites in the combined pathway - namely, 
the two pathways producing P and a pathway catabolizing P - remain decoupled. 
The simplest way to see this is to note that the process describing the synthesis of P, 
being the sum of two Poisson processes, is still a Poisson process. The pathway which 
catabolizes P is therefore statistically identical to an isolated pathway, with an incoming 
flux that is the sum of the fluxes of the two upstream pathways. More generally, the 
Poissonian nature of this process allows for different pathways to dump or take from 
common metabolite pools, without generating complex correlations among them. 

3. 5. Converging pathways - reaction with two fluctuating substrates 

As mentioned above, some reactions in a biosynthesis pathway involve side-reactants, 
which are assumed to be abundant (and hence at a constant level). Let us now discuss 
briefly a case where this approach fails. Suppose that the two products of two linear 
pathways serve as precursors for one reaction. This, for example, is the case in the 
arginine biosynthesis pathway, where L-ornithine is combined with carbamoyl-phosphate 
by ornithine-carbamoyltransferase to create citrulline (figure [3](b)) [24J. Within a flux 
balance model, the net fluxes of both substrates must be equal to achieve steady state, 
in which case the macroscopic Michaelis-Menten flux takes the form 
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Figure 4. Time course of a two-substrate enzymatic reaction, as obtained by a 
Gillespie simulation [U]. Here c = 3£ _1 , £;+ = 5i _1 and fc_ = 2t~ 1 for both substrates, 
t being an arbitrary time unit. 



Here [Si^] are the steady-state concentrations of the two substrates, and -Kjv/1,2 the 
corresponding MM-constants. However, flux balance provides only one constraint to a 
system with two degrees of freedom. 

In fact, this reaction exhibits no steady state. To see why, consider a typical time 
evolution of the two substrate pools (figure @J. Suppose that at a certain time one of 
the two substrates, say Si, is of high molecule-number compared with its equilibrium 
constant, mi ^> K\. In this case, the product synthesis rate is unaffected by the precise 
value of mi, and is given approximately by v max m 2 / (m 2 + K 2 ). Thus, the number m 2 of 
5*2 molecules can be described by the single-substrate reaction analyzed above, while mi 
performs a random walk (under the influence of a weak logarithmic potential), which is 
bound to return, after some time r, to values comparable with K\. Then, after a short 
transient, one of the two substrates will become unlimiting again, and the system will 
be back in the scenario described above, perhaps with the two substrates changing roles 
(depending on the ratio between K\ and K 2 ). 

Importantly, the probability for the time r during which one of the substrates is 
at saturating concentration scales as r~ 3 / 2 for large r. During this time the substrate 
pool may increase to the order yfr. The fact that r has no finite mean implies that this 
reaction has no steady state. Since accumulation of any substrate is most likely toxic, 
the cell must provide some other mechanism to limit these fluctuations. This may be one 
interpretation for the fact that within the arginine biosynthesis pathway, L-ornithine is 
an enhancer of carbamoyl-phosphate synthesis (dashed line in figure 0(b)). 

In contrast, a steady-state always exists if the two metabolites experience linear 
degradation, as this process prevents indefinite accumulation. However, in general one 
expects enzymatic reactions to dominate over futile degradation. In this case, equal 
in-fluxes of the two substrates result in large fluctuations, similar to the ones described 
above [31] • 



Stochastic fluctuations in metabolic pathways 



14 



4. Discussion 

In this work we have characterized stochastic fluctuations of metabolites for dominant 
simple motifs of the metabolic network in the steady state. Motivated by the analogy 
between the directed biochemical pathway and the mass transfer model or, equivalently, 
as the queueing network, we show that the intermediate metabolites in a linear pawthway 
- the key motif of the biochemical netrowk - are statistically independent. We then 
extend this result to a wide range of pathway structures. Some of the results (e.g., 
the directed linear, diverging and cyclic pathways) have been proven previously in 
other contexts. In other cases (e.g., for reversible reaction, diverging pathway or with 
leakage/dilution), the product measure is not exact. Nevertheless, based on insights 
from the exactly solvable models, we conjecture that it still describes faithfuly the 
statistics of the pathway. Using the product measure as an Ansatz, we obtained 
quantitative predictions which turned out to be in excellent agreement with the numerics 
(figure [2]). These results suggest that the product measure may be an effective starting 
point for quantitative, non-perturbative analysis of (the stochastic properties) of these 
circuit /networks. We hope this study will stimulate further analytical studies of 
the large variety of pathway topologies in metabolic networks, as well as in-depth 
mathematical analysis of the conjectured results. Moreover, it will be interesting to 
explore the applicability of the present approach to other cellular networks, in particular, 
stochasticity in protein signaling networks |2j, whose basic mathematical structure is also 
a set of interlinked Michaelis-Menton reactions. 

Our main conclusion, that the steady-state fluctuations in each metabolite depends 
only on the properties of the reactions consuming that metabolite and not on fluctuations 
in other upstream metabolites, is qualitatively different from conclusions obtained for 
gene networks in recent studies, e.g., the "noise addition rule" [Til [15] derived from the 
Independent Noise Approximation, and its extension to cases where the singnals and 
the processing units interact [IT] . The detailed analysis of [17], based on the Linear 
Noise Approximation found certain anti-correlation effects which reduced the extent of 
noise propagation from those expected by "noise addition" alone [TH [15] . While the 
specific biological systems studied in [T7J were taken from protein signaling systems, 
rather than metabolic networks, a number of systems studied there are identical in 
mathematical structure to those considered in this work. It is reassuring to find that 
reduction of noise propagation becomes complete (i.e., no noise propagation) according 
to the analysis of [17], also, for Poissonian input noise where direct comparisons can 
be made to our work (ten Wolde, private communication). The cases in which residue 
noise propagation remained in [17J, corresponded to certain "bursty" noises which is non- 
Poissonian. While bursty noise is not expected for metabolic and signaling reactions, it 
is nevertheless important to address the extent to which the main finding of this work 
is robust to the nature of stochasticity in the input and the individual reactions. The 
exact result on the cyclic pathways and the numerical result on the directed pathway 
with feedback inhibition suggest that our main conclusion on statistical independence 
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of the different nodes extends significantly beyond strict Poisson processes. Indeed, 
generalization that preserve this property include classes of transport rules and extended 
topologies [3?t 138]. 

The absence of noise propagation for a large part of the metabolic network allows 
intermediate metabolites to be shared freely by multiple reactions in multiple pathways, 
without the need of installing elaborate control mechanisms. In these systems, dynamic 
fluctuations (e.g., stochasticity in enzyme expression which occurs at a much longer 
time scale) stay local to the node, and are shielded from triggering system-level failures 
(e.g., grid- locks). Conversely, this property allows convenient implementation of controls 
on specific node of pathways, e.g., to limit the pool of a specific toxic intermediate, 
without the concern of elevating fluctuations in other nodes. We expect this to make 
the evolution of metabolic network less constrained, so that the system can modify its 
local properties nearly freely in order to adapt to environmental or cellular changes. The 
optimized pathways can then be meshed smoothly into the overall metabolic network, 
except for junctions between pathways where complex fluctuations not constrained by 
flux conservation. 

In recent years, metabolomics, i.e., global metabolite profiling, has been suggested 
as a tool to decipher the structure of the metabolic network J39JHD]. Our results suggest 
that in many cases, steady-state fluctuations do not bare information about the pathway 
structure. Rather, correlations between metabolite fluctuations may be, for example, the 
result of fluctuation of a common enzyme or coenzyme, or reflect dynamical disorder [27J. 
Indeed, a bioinformatic study found no straightforward connection between observed 
correlation and the underlying reaction network |41j. Instead, the response to external 
perturbation [281 139| [42] may be much more effective in shedding light on the underlying 
structure of the network, and may be used to study the morphing of the network under 
different conditions. It is important to note that all results described here are applicable 
only to systems in the steady state; transient responses such as the establishment of 
the steady state and the response to external perturbations will likely exhibit complex 
temporal as well as spatial correlations. Nevertheless, it is possible that some aspects of 
the response function may be attainable from the steady-state fluctuations through non- 
trivial fluctuation-dissipation relations as was shown for other related nonequilibrium 
systems [22] 133]. 
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Supporting Material 

Appendix A. Microscopic model 

Under the assumption of fast equlibration between the substrate and the enzyme, the 
probability of having N$e complexes given m substrate molecules and Ne enzymes is 
given by equation (3) of the main text. To write the partition function explicitly, we 
define u(x) = U{x, 1 — m — Ne', —K), where U denotes the Confluent Hypergeometric 
function [15] . One can then write the partition sum as Z m<NE = (—K)~ NE u(—m). 
The turnover rate is then given by w m = y E [—mu(l — m)]/[u(— m)], which can be 
approximated by Equation (4). 



Appendix B. Influx of metabolites 



A metabolic reaction in vivo can be described as turnover of an incoming flux of substrate 
molecules, characterized by a Possion process with rate c, into an outgoing flux. To find 
the probability of having m substrate molecules we write down the Master equation, 



d 

—7i(m 
dt 



) = [c(a — 1) + (a — l)w m ] Tc(m) = c[7r(m—l)—7r(m)]+[w m+1 7r(m+l)—w m n(m)] , (B.l 



where we took the opportunity to define the lowering and raising operators a and a, 
which - for any function h(n) - satisfy ah(n) = h(n—l), ah(0) = 0, and ah(n) = h(n+l). 
The first term in this equation is the influx, and the second is the biochemical reaction. 
The solution of this steady state equation is of the form 7r(m) ~ c m / [Xb=i Wk (up to a 
normalization constant), as can be verified by plugging it into the equation, 

1) .\ fn(m+l) 



TT(m 



7r(m) 



7r(m) 



-w m +\ — w m 



Wm 
C 



- 1 



w 



W m+ i - w, r 



m+1 



Using the approximate form of w m , as given in (4), the probability 7r(m) takes the 
form, 



7r m 



m + K+(N E -l] 



m 



(1-z) 



K+N Ez m 



(B.3) 



as given in equation (5) of the main text. 



Appendix C. Directed linear pathway 



We now derive our key results, equation (8) (The result has been derived previously in 
the context of queueing networks [IB] , and of mass-transport systems [17] ) . To this end 
we write the Master equation for the joint probability function tc = 7r(m.i, m.2, • • • , tul), 

L-l 



d_ 

dt 



71 



c(ai - 1) + ^(dia i+ i - + (« L _ i) w 



IIIL 



71 



(C.l) 



which generalizes ( IB. II ). As above, a, and &i are lowering and raising operators, acting 
on the number of Si molecules. The first term in this equation is the incoming flux c 
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of the substrate, and the last term is the flux of end product. Let us try to solve the 
steady-state equation by plugging a solution of the form 7r(mi, m 2 , ■ • • , m L ) = F] 9i{Tni), 
yielding 

ftCmi-l) ^.yvfi) giK + l)^+i(m i+ i - 1) (i) (L) ff L (m L + 
9i\ m i) ~[ gi{mi)gi+i{mi+i) 9l\T^l) 

Motivated by the solution to ( IB. II ). we try to satisfy this equation by choosing 
<7i(m) = c m j YYk=i w k • With this choice we have g(m + l)/g(m) = c/w m+ i and 
g(m — l)/g(m) = w m /c. It is now straightforward to verify that indeed 



w 



(1) \ L-1 



mi 



- 1 + e («sa - -s) + ( c - o = ° • ( c - 3 ) 
/ i=i 

Finally, in our choice of gi(m) we replace by the MM- rate ViUii/ '(jrii + Ki), and find 
that in fact gi(m) = 7Tj(m), namely 

L 

7r(mi,m 2 , . . . ,m L ) = J^7Ti(mi) , (C.4) 

i=l 

as stated in (8). 



Appendix D. End-product inhibition 

Equation (13) of the main text is a self-consistent equation for the steady- state flux c 
through a pathway regulated via end-product inhibition. Using considerations analogous 
to what led to the exact result on the product measure distribution for the cyclic 
pathways, we conjecture that even for the present case of end-product inhibition, the 
distribution function can still be approximated by the product measure ( IC.4I ) with the 
form of the single node distributions given by ( IB. 3 1 ). The flux c enters the calculation of 
the average on the right-hand side through the probability function 7r(m). Solving this 
equation for c yields the steady state current, and consequently determines the mean 
occupancy and standard deviation of all intermediates. 

To verify the validity of this conjecture, and to demonstrate its application, we 
consider the case h = 1. In this case one can carry the sum, and find 

oo 

c= J2 c [l + (m L /Ki) h ]~\ L {m L ) (D.l) 

m L =0 

= c (l - z) K ^F x {Kj, K L] Kj + 1; z) 

with z = c/vl and 2F1 the hypergeometric function [45]. This equation was solved 
numerically, and plotted in supporting figure lD2l a) for some values of Ki and Kl. Note 
that predictions based on the product measure (lines) are in excellent agreement with 
the results of numerical simulation (circles) for the different sets of parameters tried. 
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Figure Dl. The steady-state distribution 7r(m) of a metabolite, that experiences 
enzymatic reaction (with rate w m — vm/(m + K — 1)) and linear degradation (with 
rate /3to), as given by equation (10) of the main text. Here K — 100 and v = 2cq. 



Results obtained from equation ( ID. II ) can be used, for example, to compare the 
flux that flows through the noisy pathway with the mean- field flux cmf, obtained when 
one ignores fluctuations inmi, i-e., 

CMF = 1 + (s°/K,) h ' (D ' 2) 
The fractional difference 5c = (c — cmf)/cmf is plotted in supporting figure lD2l (b). The 
results show that number fluctuations in the end-product always increase the flux in the 
pathway since 5c > always. Quantitatively, this increase can easily be several percent. 
For large Cq, a simplifying expression can be derived by using an asymptotic expansion 
of the hypergeometric function [45J. For example, when Kj < Kl, 

(1 - zf^F^K;, K L - Kj + 1; z) ~ - V ^ Kl , (D.3) 

1 + K L — Ki 

which yields 

C ~ CMF „ J- V Jl . (D.4) 
cmf Ki c 

Thus the effect of end-product fluctuations on the current is enhanced by stronger 
binding of the inhibitor (smaller Ki), as one would expect. We note that obtaining 
these predictions from Monte-Carlo simulation is rather difficult, given the fact that 
one is interested here in sub-leading quantities. 
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Figure D2. Pathway with end-product inhibition. The influx rate is taken to be 
co/(l + rriL/Kj), and thus the steady-state flux is given by equation (12) of the main 
text, with h = 1. (a) Assuming that different metabolites in the pathway remain 
decoupled even in the presence of feedback regulation, (12) can be approximated 
by ( ID. II ). Numerical solutions of equation ( ID. II ) (lines) are compared with Monte- 
Carlo simulations (symbols). Values of parameters are chosen randomly such that 
100 < Ki < 1000 and c < < 10c. For the data presented here, vl = 2.4c. We 
find that ( ID. II ) yields excellent prediction for the steady-state flux, (b) Neglecting 
fluctuations altogether yields a mean- field approximation for the flux, cmf, given in ( 
ID. 21 ). For the same data of (a), we plot the fractional difference 5 C — (c — cmf)/c We 
find that steady-state flux is increased by fluctuations, and thus taking fluctuations 
into account (even in an approximate manner) better predicts the steady-state flux. 
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